First-order phase transitions in two-dimensional off-lattice liquid crystals 
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We consider an off-lattice liquid crystal pair potential in strictly two dimensions. The potential 
is purely repulsive and short-ranged. Nevertheless, by means of a single parameter in the potential, 
the system is shown to undergo a first-order phase transition. The transition is studied using mean- 
field density functional theory, and shown to be of the isotropic-to-nematic kind. In addition, the 
theory predicts a large density gap between the two coexisting phases. The first-order nature of 
the transition is confirmed using computer simulation and finite-size scaling. Also presented is an 
analysis of the interface between the coexisting domains, including estimates of the line tension, as 
well as an investigation of anchoring effects. 

PACS numbers: 64.70.Md, 02.70.-c, 61.30.Cz, 61.30.Hn 



I. INTRODUCTION 

The phase behavior of liquid crystals in two dimen- 
sions continues to be an interesting topic. On the one 
hand, at least for lattice liquid crystals, there is a clear 
resemblance to the planar spin or XY-model In 
fact, the Lebwohl-Lasher (LL) model 0, which is one 
of the standard liquid crystal models, maps exactly onto 
the XY-model in two dimensions. The XY-model does 
not support long-range order [sj and, consequently, long- 
range nematic order is believed to be absent in many 
two-dimensional (2D)liquid crystals as well [1, @|, Q (an 
exception being Ref. |7i) . In addition, the XY-model fea- 
tures a Kosterlitz-Thouless (KT) transition Conse- 
quently, phase transitions in two-dimensional liquid crys- 
tals are often interpreted in terms of the KT scenario 
0) B B [13) [HI- In particular, the KT transition is a 
continuous transition, as opposed to first-order. As a re- 
sult, the possibility of a first-order transition occurring in 
a two-dimensional liquid crystal, has received relatively 
little attention. 

Interestingly, computer simulations of appropriately 
generalized XY-models have shown that the possibility 
of also a first-order transition occur ring in these systems 
should be taken seriously [H, [H, Il4| More recently, 
these findings were put on firm mathematical ground by 
van Enter and Shlosman [l^ [1^ [13 • In particular, it 
was demonstrated that Hamiltonians of the form 
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+ LUi ■ Lb. 



(1) 



undergo first-order phase transitions when p is large [17| . 
In Eq. ([T]) , the sum is over nearest neighbors on a square 
lattice, and uji are two-dimensional unit vectors. The 
usual XY-model (up to a trivial constant and multiplica- 
tive factor) is recovered for p = 1; the generalization is to 
also consider p ^ 1. In computer simulations, the first- 
order transition was already observed at p = 50 [lH . Of 
course, for small p, the KT-scenario is ultimately recov- 
ered again, for which the transition is continuous. 



The results obtained for generalized XY-models have 
similar consequences for two-dimensional liquid crystals. 
This was recently demonstrated in Ref. [H using a gen- 
eralized version of the LL-modcl 



Ti-gLL — — 



E 



(2) 



the essential difference with Eq. ([T]) being inversion sym- 
metry under uji <-» —ilii of the particle orientations. In 
particular, it was shown that Eq. ^ undergoes a first- 
order temperature-driven transition, from an isotropic to 
a quasi-nematic phase, provided p is sufficiently large. 
Again, the threshold value is at p w 50 In the 

isotropic phase, the orientational correlations decay ex- 
ponentially; in the quasi-ncmatic phase, they decay alge- 
braically. Both phases thus lack long-range order in the 
thermodynamic limit, in agreement with the Mermin- 
Wagner theorem 0]. Consequently, the nematic order 
parameter cannot be used to describe the first-order tran- 
sition in Eq. ([2]). Instead, a valid order parameter is 
the energy density, which shows a "jump" at the tran- 
sition temperature. By keeping the energy density fixed 
at some value in the coexistence region, phase coexis- 
tence between isotropic and quasi-nematic domains can 
be realized [18]. 

The aim of this paper is to investigate how robust these 
findings are when also off-lattice liquid crystals in two di- 
mensions are considered. To this end, we formulate an 
off-lattice model of a liquid crystal, which is somewhat 
inspired by the lattice Hamiltonian of Eq. ^ . The model 
will be presented in Section [Til Next, the phase behav- 
ior of this model is studied, using theory and simulation. 
Indeed, both the theory and the simulation find strong 
evidence for the existence of a first-order transition, in- 
cluding a pronounced coexistence region. The coexis- 
tence region will be analyzed in some detail, including 
estimates of the line tension between the coexisting do- 
mains. An important improvement over the lattice model 
of Eq. ([2]) is that the transition in the off-lattice model 
is characterized by a "jump" in the particle density. In 
other words, phase coexistence can now be studied by 



2 



keeping the overall particle density fixed at some appro- 
priate value, rather than the energy density. This finding 
is relevant for possible experiments, where the condition 
of fixed density is rather easy to implement (in contrast, 
keeping the overall energy fixed in an experiment would 
be much more difficult). At the same time, we find that 
the analogue of the p-exponent in Eq. ^ must be quite 
large, before the first-order transition begins to show-up. 
Whether such "sharp" interactions can be realized exper- 
imentally is not yet clear, but some suggestions are made 
toward the end of this paper. 



II. OFF-LATTICE LIQUID CRYSTAL IN TWO 
DIMENSIONS 



In this paper, we consider an ensemble of particles, 
whose positions and orientations are confined to a two- 
dimensional plane. The particles interact with each other 
via a pair potential Vij of the form 



(1 



\UJt ■ LUj\ 



u(r). 



1 



(3) 
(4) 



with r.y = Vj — n, r — \rij\, —1/2 < v < 1/2, and cou- 
pling parameter e > 0. In what follows, factors of k^T 
are absorbed into the parameter e, with T the tempera- 
ture, and fee the Boltzmann constant. The quantity uJi 
is a two-dimensional unit vector denoting the orientation 
of the i-th particle; is a two-dimensional vector de- 
noting the coordinate of the center of mass of the i-th 
particle. The radial function u{r) in Eq. is assumed 
to be strictly positive and short-ranged. Here, we take a 
simple step function 



u{r) = 



1 r < a, 
otherwise. 



(5) 



with a the particle diameter which will henceforth serve 
as our unit of length. 

As for the lattice Hamiltonian of Eq. 0, the poten- 
tial is constructed such that inversion symmetry is main- 
tained. Note also that Vij is purely repulsive. Neverthe- 
less, we expect a first-order phase transition to take place, 
either at low temperature, or at high density, provided p 
is large. The purpose of the parameter ly is to introduce 
a coupling between the orientational and translational 
degrees of freedom. By setting v = 0, no such coupling 
occurs, and the potential becomes separable. For this 
special case, Straley has rigorously proved the absence 
of long-range nematic order (19j . As far as we know, the 
absence or presence of long-range order for the case ^ 
is still an open question. 



III. DENSITY FUNCTIONAL THEORY 

Within density functional theory (DFT) the thermo- 
dynamics and structure of a fluid is described by a func- 



tional il[p] of the one-particle distribution p{r,uj). The 
density functional is such that it is minimized for a given 
{fi, A, T) by the equilibrium one-particle distribution and 
the minimum value of the functional is the grand poten- 
tial [20] (in the present study in two dimensions, A is the 
system area) . Here we use a simple mean-field functional 
which, in the absence of an external potential, can be 
cast in the following form 



n[p\= J drda)p(r, lo) (ln[Vp(r, w)] - 1 - fi) (6) 
+ i J drd(jj J dr' doj' v{Ar\uj, Lu')p{r, uj)p{r' ,u}'), 

with expressed in units fceT, Ar = r' — r, V the 2D 
thermal volume of the particle and p the chemical po- 
tential. The functional is known to give accurate results 
for dense fluids of soft spheres with bounded potentials 
[2ll . [23 , [2^ . Recently, it was shown that the approach 
also works well for fluids of soft anisometric particles [13] . 
The minimum condition on the functional, SQ[p\/Sp = 0, 
leads to a nonlinear integral equation 

ln[Vp(r, u})] + [ dr'dLj'v{Ar; Cj, uo')p{r' , Cj') = p, (7) 



to be solved for the equilibrium distribution p(r',cD') at 
a given p. 



A. Bulk phase diagram 

Let us first focus on localizing the isotropic-to-nematic 
transition for bulk systems. As the average density in 
both states is spatially uniform, i.e. independent of r, 
we may write p(r,cD) = pf{uj), where p is the bulk den- 
sity and / an orientational distribution, subject to the 
normalization condition J dcuf^Lu) = 1. In the isotropic 
(I) state, all orientations are equally probable and /i is a 
constant while in the nematic (N) phase /n is expected 
to be strongly peaked around the nematic director which 
we assume to be spatially uniform. This implies that 
long-range orientational order is always present in theory 
since fluctuations or local defects in the director field are 
not taken into account. Introducing the angle € [0,7r] 
between the particle and the nematic director we may 
rewrite Eq. ^ into a self-consistency equation for f{ip) 



exp 



p / d^'E{^,^')fy) 



(8) 



where Z = J d(y9exp[- • ■] to ensure normalization. Note 
that /(</?) — /(tt — Lp) due to the inversion symmetry. 
The kernel £^(1^5, <^') is given by the following spatial in- 
tegration 



E{if,ip') = j dArt>(Ar 



UJ,Uj'). 



(9) 
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The spatial integration in Eq. ([9]) can then be carried out 
without difficulty to give 



E{lp, If') = e(l + v)T:a^ (1 - | cos((p' - ip)\P) . 

With this result it is expedient to rewrite Eq. ([ 
following way 



(10) 



in the 



/(^) =Z-iexp 



c / V|cos(^'-'^)r/(^') 



(11) 



where we have introduced the effective dimensionless den- 
sity 



7r(l + v)pa^e. 



(12) 



The combination pe in Eq. \12\ illustrates the fact that 
phase transitions in our model can be brought about by 
either increasing the density or lowering the temperature. 
It is clear that the isotropic distribution // = I/tt is a 
trivial solution of Eq. pT|) at all densities c. However 
at high densities non-trivial, i.e. nematic solutions are 
expected to branch off from the isotropic one. To locate 
the branching point (denoted by c*) we will perform a 
simple stability analysis of Eq. pT|) . along the lines of 
Refs. m, [2^. Let us introduce the following expansions 
around the isotropic solution 



f{ip) — — [1 + aai cos2(^ -I- a^a2 cos4(^ -t- • ■ •] 

TT 

c = Co + aci + a^C2 + • • ■ , 



(13) 



in terms of a single order parameter a. Likewise, we may 
expand the kernel as follows 

cos{ip' — 1^9)1^ — k2n{p) cos{2n(p) cos(2n(^'), (14) 



with coefficients 



k2n{p) 



(15) 



— I d(p I d(p'\coa{(p' — (p)\P coa{2n(p) cos{2rnp'). 



Inserting all expansions back into Eq. (fTTI) and keeping 
all contributions up to 0{a) gives the branching or bi- 
furcation point 



Co = 2/fc2(p), 



(16) 



implying that stable nematic solutions are expected for 
c > c*. We remark that for p = I the branching density 
c* = 37r/2 is identical to that of the 2D Onsager theory 
for infinitely thin hard needles [l^, H^l • 

To verify the thermodynamic stability of the nematic 
solutions close to the branching density we must analyze 
its free energy. The dimensionless Helmholtz free energy 
(ignoring all irrelevant constants) is given by 

F[f]/N ^ (ln[c^/(^)]) + ^(((1 - I cos(^' - ^)r))), (17) 



where (• • ■ ) = / d(pf{(p). The free energy difference 
AF = Fn — Fi between the nematic and isotropic states 
close to the branching point can be written in Landau 
form in terms of the nematic order parameter a 



AF/N 



(18) 



The coefficients can be obtained from extending the bi- 
furcation analysis to higher order in a [using the expan- 
sions Eq. (fT3|) and Eq. (fT4|) ] and performing an order- 
by-order solution of Eq. (fTTj) up to 0{a'^). The algebra 
is straightforward but tedious and we will only give the 
final outcome here. It turns out that A,B ~ Q while 



C=- 



k2{p)-2h{p) 



(19) 



The order of the transition depends on the sign of C. 
If it is negative, small nematic perturbations around the 
branching point immediately stabilize the system and the 
transition is continuous. If C is positive the incipient 
free energy difference goes up for small a which means 
that the actual phase transition must involve a density 
jump and is first-order. In the latter case the coexist- 
ing densities are indicated by binodal curves which can 
be computed in the usual way by requiring the pres- 
sure and chemical potential to be equal in the coexist- 
ing phases. The equilibrium f{ip) for a given nematic 
density c is obtained numerically from the consistency 
equation Eq. pTjl by dividing the interval [0,7r/2] into 
100 equidistant grid points and employing the iteration 
scheme outlined in Ref. [H. 

Fig. [1] shows that the isotropic-nematic transition is 
continuous for small p but becomes first-order for p > 8. 
The full phase diagram in Fig. [2] shows that the actual 
crossover from continuous to first-order (marked by the 
tri-critical point where all phase lines meet) is located at 
a somewhat lower value for p, namely p — 4.7. The dis- 
crepancy is expected since the bifurcation analysis usu- 
ally provides an upper estimate for the critical points. 
We also remark that the non-monotonic behavior of the 
critical density at low p is consistent with the simulation 
results of the generalized LL-model reported in Ref. [H. 
The steep increase of the nematic order parameter in 
Fig. [3] suggests that a considerable degree of nematic or- 
der is expected in the coexisting nematic phase at large p. 

1. Asymptotic results for large p 

For very large p (say larger than 100) the solution of 
Eq. (jlip on a grid becomes numerically awkward since 
f{(p) gets extremely peaked at </5 = (and tt). It is 
therefore tempting to formulate a simple variational the- 
ory that allows us to access the phase diagram at asymp- 
totically large p. Indeed, for highly nematic states f{ip) 
is well-described by a gaussian trial function with varia- 
tional parameter /3 ^ 1 p^ . 

fm ^ [ — ] exp 



for < »> < 2 . (20) 
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and its mirrored version f{TT~Lp) for the interval tt/2 < 
ip < TT. Inserting the gaussian into Eq. pT|) and integrat- 
ing yields the following asymptotic result for the ideal 
free energy (first term) 



(In[c7r/(v5)]) 



In c 
In c, 



1 



■ In 2tt(3 



1 

2' 



(N) 
(I) 



valid for large /3. Using the approximation 

cosxl^' w exp[-(l/2)pa;^], (p > 1) 



(21) 



(22) 



in the excess free energy (second term) allows us to cal- 
culate the double orientational averages analytically 



((icos(^'-^)n) 



f3 



f3 + 2p 



1/2 



(N) 

(I) (23) 



The value of /3 is fixed (at a given density) by minimizing 
the total nematic free energy with respect to the varia- 
tional parameter. Some rearranging then leads to the 
following minimization condition 



2) -c'^ = 0, 



(24) 



which has to be solved numerically along with the coex- 
istence equations for the chemical potential and pressure 
at a given p. These follow straightforwardly from the 
free energy Eq. pT]) . The resulting binodals and coexis- 
tence chemical potential are shown in Fig. [3] and Fig. O 
respectively. 

Another feature we want to point out is that the bi- 
furcation density c* at large p scales as c* oc p^^^. The 
scaling relation can be easily established by making an 
asymptotic expansion of ^2 from Eq. (jlSp . For the regime 
p ^ 1 one can show with the aid of Eq. that 



k2ip) K / d((^'-^)exp[-(l/2)p((^'-(^)2 
Jo 



(25) 



up to leading order, hence c* oc p^/^. This result is anal- 
ogous to the scaling of the critical coupling constant in 
the LL-simulations of Ref. llSl. 



B. Isotropic- nematic interface 

To assess the properties of the interface between the 
coexisting isotropic and nematic phases we have to go 
back to our initial DFT formulation in Section IIIII If we 
assume the interface to be flat with a surface normal x, 
the one-particle density will be non-uniform along this 
direction and depend on the spatial coordinate a; = r • a; 



o 




FIG. 1: Landau coefficient C [Eq. (fT9|) ] versus p. For p > 8 
the isotropic-nematic transition is first-order; for smaller p, it 
is continuous. 




FIG. 2: Theoretical phase diagram. The dashed curve repre- 
sents the nematic branching line c* , calculated from Eq. (|16|l . 
The binodals are given by the solid curves. A tri-critical point 
is located at p = 4.7. 



and angle tp. It is a solution of the integral equation 
Eq. dll) with H = H*, 

ln[Vp(x,vp)]+ (26) 
dx'dp'E,{Ax; if, p';d)p{x', p') = 



fi* being the chemical potential at coexistence, and Aa; = 
x' — X, subject to the boundary conditions 

p{x,(p) = pifi, {x ~oo), 

p{x,p) = pNfNip), (x^oo), (27) 



(recall that p = arccos(n ■ w)). In Eq. P^ . the kernel 
Ex is defined as the pair potential (at fixed orientations) 
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FIG. 3: Nematic order parameter S = (cos 2ip) corresponding 
to the nematic binodal in Fig. [S] 
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FIG. 5: Coexistence chemical potential /i* versus p corre- 
sponding to Fig. |4] The solid curve is the theoretical result; 
squares are simulation results. 
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FIG. 4: Phase diagram of the two-dimensional liquid crys- 
tal model of Eq. ([3} using e = 2.5. The solid curves show 
the theoretical binodals, obtained using the large p approxi- 
mation. Plotted is the dimensionless density p = Nc? jA of 
the isotropic phase (left curve) and the nematic phase (right) 
versus p. Also shown are the corresponding simulation results 
(squares), where the dashed lines serve to guide the eye. 



averaged over the distance Ay perpendicular to i, along 
which the system is homogeneous 



/oo 
-OO 



(28) 



Note that because of the broken spatial symmetry the 
integral depends implicitly on the anchoring angle -d = 
arccos(ri • x) between the nematic director h and the sur- 
face normal. This becomes manifest when we focus on the 
translation-rotation coupling contribution a in Eq. ([?]). 
In explicit form it reads 



with Af = {Ax, Ay}/(Aa; 
difference unit vector, a) = 
tation matrix 



' + Ay^)-'^/^ the center-of-mass 
{cos iy9, sin 93} and TZ^ the ro- 



cos -d ~ sin -d 
sin -d cos "d 



(30) 



Clearly, if = there is no dependence on d and the 
interfacial profiles are identical for all anchoring angles 
30]. \i V ^ Q the translational and rotational degrees 
of freedom are coupled and the interfacial properties will 
in general be dependent upon the anchoring angle. In 
particular we are interested in the line tension 7 which 
can be extracted from the equilibrium interfacial density 
profile. Inserting the one-body density into the functional 
Eq. ([7]) yields the minimum grand potential 



^min — 

dx I dipp{r,tp) 



(31) 



2MVp(r,^)] 



1 - -a* 
2^ 



The line tension 7 [/i* (p) , i?] is then obtained from the 
standard thermodynamic relation 7 — (rJmin + PA)/L, 
with P the coexistence pressure. 

In Fig. [5] we show the line tension for 1^ = 0. To facili- 
tate comparison with simulations later on we will hence- 
forth fix the coupling parameter to e = 2.5. Note that, 
owing to Eq. (fT2| , changing this value does not give qual- 
itatively different results but merely constitutes a linear 
shift in 7 and p. The increase of 7 as a function of p re- 
flects the transition becoming strongly first-order at large 
p. The corresponding interfacial profiles for the density 
p{x) and the nematic order parameter S{x), defined as 



a = 1 + 1^ [(Af ■ 7^^cZ))2 + (Af • 7^^>a)')^] 



(29) 



p{x) = J dipp{x,(p), 

S{x) — p^^{x) / dipp{x,tp) cos2ip, (32) 
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FIG. 6: Line tension 7 (in units k^T /a) for e = 2.5 versus p. 
(a) DFT result, (b) Computer simulation results (squares), 
together with the DFT result (solid curve) for comparison. 
Note the different p-range in the above plots. 
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FIG. 7: Interfacial profiles for (a) the density p(x) and 
the nematic order parameter S{x) for 1^ = 0. 



(b) 



are shown in Fig. [71 As expected, the interface becomes 
sharper for large p. For p — 100 small density oscillations 
on the isotropic side occur, which point to a layering 
effect induced by the interface. 

Let us now focus on the anchoring behavior for a fixed 
value of p = 75 and v ^ 0. The results, summarized in 
Fig. [HI show a strong dependency of the line tension on 
the anchoring angle t?, especially for j/ < 0. Recall that 
end-to-end pair configurations are energetically more fa- 
vorable than side-to-side ones in this case. While for 

V < Q the minimal tension is at ■)? = implying perpen- 
dicular surface anchoring, a positive v seems to be asso- 
ciated with parallel anchoring. Investigations for other p 
reveal that this phenomenon is robust. Microscopic in- 
formation extracted from the profiles for v — —0.45 is 
shown in Fig. [SI These reveal that the layering effect is 
considerably influenced by the anchoring angle. In par- 
ticular, particles show enhanced localization across the 
interface if the nematic director is forced to be parallel 
to the interface. 

An even more dramatic effect is encountered if we lower 

V —0.5, see Fig. (TO] Note that, for v = —0.5, the end- 
to-end configurations have zero repulsion and are there- 
fore strongly favored. This explains the tendency of the 
system to form string like clusters along the interface as 
indicated by the sharp density peaks in Fig.[Tni However 
v = —0.5 seems to be a rather pathological case. The 
density modulations penetrate deeply into the isotropic 
bulk which raises strong suspicions as to whether the 
isotropic fluid state is stable against clustering or crystal- 
lization. It remains to be checked by simulations whether 
the IN transition is really pre-empted by a freezing tran- 
sition in this case. 
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FIG. 8: Line tension 7 (in units k-aT/a) versus anchoring 
angle 1? (in degrees), for several values of v as indicated. 



IV. COMPUTER SIMULATION 



Next, we will confirm some of the theoretical results 
using computer simulation. To this end, grand canonical 
(GC) Monte Carlo simulations fsH, Is^l are performed, at 
chemical potential /i, for the model of Eq. ([3]) using e — 
2.5. Again, for the radial part, we take the step function 
of Eq. ([5]). In this work, we restrict the simulations to the 
case J/ = in Eq. ([3]). A computer simulation study of 
anchoring effects, predicted by the theory to occur when 
J' 7^ 0, will be postponed to a future publication. 
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FIG. 9: Interfacial profiles for p = 75, v = —0.45 at difTerent 
anchoring angles j? (in degrees). 
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FIG. 10: Interfacial profiles for p = 75 and anchoring angle 
1? = 90 degrees. 



A. Finite-size scaling results for p = 75 

We begin the simulations using p = 75 in Eq. ([3]). 
Inspired by our theoretical results, we expect that, for 
p = 75, a first-order phase transition should occur, when 
the chemical potential is tuned to its coexistence value 
fi* . In GC simulations, /i* is determined from the distri- 
bution Pl fi) , defined as the probability to observe a 
particle density p in the system. In order to also simulate 
the regions of low probability, a biased sampling scheme 
is implemented (ssl . [33 | . The distribution {p, p) de- 




FIG. 11: (a) Coexistence distributions Pl{p, fi) for various 
system sizes L as indicated, with p = Na^/A the dimension- 
less particle density. The coexistence between the two phases 
is manifested by "equal-area" under the peaks, (b) The corre- 
sponding distributions Pl (S) of the nematic order parameter 
S, given by Eq. (fM)l . 



pends on the chemical potential p, as well as on the size 
of the L X L simulation square (periodic boundary con- 
ditions are assumed). At p — p* , Pl{p,p) becomes bi- 
modal, with two peaks of equal area. In computer simu- 
lations, the transition can thus be located by varying p 
until the "equal-area" criterion is obeyed. 

Indeed, we find that bimodal density distributions can 
be realized in this way. Fig. Illf a) gives some examples 
for various system sizes L. By increasing L, the peaks be- 
come more narrow. This is to be expected because, in the 
thermodynamic limit L — > oo, one has a sharp transition, 
and, consequently, a distribution featuring two (5-peaks. 
In addition, closer inspection of Fig. [TlT a) also reveals a 
finite-size effect in the peak positions. Finite-size effects 
at first-order phase transitions have received consider- 
able attention [sl, [sl, [13] . We believe that the current 
state-of-the-art in this field is the rigorous treatment of 
Borgs and Kotecky [s^l- More precisely, for two-phase 
coexistence data obtained using the "equal-area" rule, 
an exponential L-dependence is predicted 



i?=|Xi-Xoo|<0(e-"^) 



(33) 



Here, is the property of interest as obtained in the 
finite system of size L, X^o the corresponding value in 
the thermodynamic limit, and r a constant. 

In Fig. [T2y a). we show the position p^ of the high- 
density (nematic) peak of Pl(p, p) as a function of L. 
Next, we use Eq. (|33)) to estimate the density poo of the 
nematic phase in the thermodynamic limit. To this end 
we have plotted, in Fig. fl^ b). the logarithm of i? as a 
function of system size L. Here, p^o was taken to be a 
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FIG. 12: Finite-size scaling ol the position ol the high-density 
peak in Pl{p,ij.) of Fig. Illf a), (a) Peak position pL as a 
function of L. The horizontal line marks the density poo in the 
thermodynamic limit, (b) Finite-size scaling analysis using 
Eq. (|33|l. The straight line is a linear fit; see details in text. 
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FIG. 13: Finite-size scaling of the coexistence chemical po- 
tential, (a) jj,^ as a function of L. The horizontal line marks 
the coexistence chemical potential /ij^ in the thermodynamic 
limit, (b) Finite-size scaling analysis using Eq. p3|) . The 
straight line is a linear fit; see details in text. 



fit parameter, and tuned until the graph of Ini? versus L 
became hnear. For p^o « 2.985, we find that InR indeed 
becomes hnear in L, which thus serves as our best esti- 
mate of the nematic density in the thermodynamic hmit. 
For completeness, this estimate has also been marked 
in Fig. [T^ a) as the horizontal hne. The density of the 
low-density (isotropic) peak can be obtained analogously. 



The corresponding scaling plot is qualitatively similar to 
Fig. [T2] and therefore not shown here. For the density 
of the isotropic phase in the thermodynamic limit, we 
obtain ~ 1.850. 

In addition, we also observe an L-dependence in the 
coexistence chemical potential /i^ . Recall that fi^ corre- 
sponds to the chemical potential at which "equal-area" 
in Pl{p,p) of the finite system is observed. Shown in 
Fig- fT^ a) is versus L. Again, Borgs and Kotecky 
predict an exponential L-dependence given by Eq. (|33p . 
We may therefore estimate as before, by plotting In R 
versus L, using for /ij^ that value at which the data be- 
come linear. The result is shown in Fig.fTWb'). Again, the 
data follow the straight line quite well, and we conclude 
p*^ ~ 5.6536. 

It is also of interest to consider the distribution -Pl(5) 
at coexistence, with S the nematic order parameter. 
Here, we define S as the maximum eigenvalue of the ori- 
entational tensor 



1 



N 



QafJ — ^ {'i'diadifj — Safs) : 



(34) 



with dia the a component (a = x, y) of the orientation 
Cbi of molecule i, 5af3 the Kronecker delta, N the number 
of particles, and system area A = . For an isotropic 
phase, in the thermodynamic limit, S equals zero. For 
a "true" nematic phase, i.e. with long-range order, one 
has 5 > in the thermodynamic limit. Note that we 
have normalized Eq. (j34p with the system area A, and 
not with N as is usually done. The reason to normalize 
with respect to A is that, in the GC ensemble, iV is a fluc- 
tuating quantity. Note also that, for a perfectly aligned 
phase, S becomes identical to the particle density N/A. 

Several distributions Pl{S) are shown in Fig. [TlTb'). 
We emphasize that all these distributions were obtained 
using that value of p at which "equal-area" in (p, p) 
was obtained. The data of Fig. [TlT b) are quite stun- 
ning. On the one hand, we observe a pronounced shift 
of the lower-peak toward S" — > 0. This is to be expected, 
since this peak corresponds to the isotropic phase. Note 
that the L-dependence of the isotropic peak does not 
reveal any physical information: an ideal gas of rods 
would reveal the same effect. The L-dependence of S 
in the isotropic phase is purely a numerical artifact. It 
stems from the fact that S is positive, since always the 
maximum eigenvalue of the orientational tensor is taken. 
Consequently, the distribution of S in the isotropic phase 
cannot be gaussian around S" = 0. Instead, the distribu- 
tion is "skewed" , with the peak always being located at 
S > Q. As the system size is increased, the isotropic peak 
shifts to zero, see the discussion by Eppenga and Frenkel 
[11], precisely what we observe. 

In contrast, on the scale of Fig. [TlTb). the nematic 
peak appears to be rather insensitive to the system size 
L. In fact, finite-size scaling of the nematic peak position 
suggests that a finite value Soo in the thermodynamic 
limit is maintained, see Fig. [TH Shown in Fig. [TTTa) is 
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FIG. 14: Finite-size scaling of the position of the nematic peak 
in Pl(S) of Fig. Illf bl. (a) Peak position Sl as a function of 
L. The horizontal line marks the nematic order parameter Soo 
in the thermodynamic limit, (b) Finite-size scaling analysis 
using Eq. (|33|) . The straight line is a linear fit; see details in 
text. 



the nematic peak position Sl as a function of L. The 
data are quite interesting, because they show an increase 
of nematic order with increasing system size L. Shown in 
Fig. [rirb) is the result of the finite-size scaling analysis 
using Eq. (j33p . Again, Soc was obtained by tuning, until 
the best collapse of the data onto a straight line occurred. 
For the nematic order parameter in the thermodynamic 
limit, we obtain Soo ~ 2.923. Note that, for v = 0, 
the liquid crystal potential of Eq. ^ is separable. For 
this special case, Straley has proved the absence of long- 
range nematic order in the thermodynamic limit [l9( . In 
other words, Soo should become zero, while our finite-size 
analysis, in contrast, suggests that 6*00 remains finite. For 
the XY-model in two dimensions, it has been shown that 
the decay of magnetic order with system size is so slow, 
one would need a sample "the size of Texas" [so] to see it. 
The most likely explanation is therefore that something 
similar also takes place in our liquid crystal model, and 
that the data of Fig. [ID^a) will eventually "turn-over" 
and decay to zero. 

Next, we consider the logarithm of the distributions, 
which essentially reflect minus the free energy of the sys- 
tem. Shown in Fig. [12] is In Pl{S) for several system 
sizes L. Again, we emphasize that the distributions were 
obtained using that value of fi at which "equal-area" in 
Pl(p, ^) was obtained. For each distribution In Pi(S'), 
we may read-off the average peak height AF, measured 
with respect to the minimum between the peaks. By 
increasing the size of the system, AF increases as well. 
As was shown by Binder, AF corresponds to the free- 
energy cost of having two interfaces in the system [ioj . 
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FIG. 15: Distributions In Pl{S) for various system sizes L as 
indicated. Also marked is AF for the L = 30 system. 
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FIG. 16: (a) Free energy barrier AF, extracted from the dis- 
tributions In Pl (S) of Fig. 1151 as a function of L. The straight 
line is a fit to Eq. (|35p . (b) Same as (a), but with AF ex- 
tracted from \n Pl{p, fJ.)- 



More precisely, in two dimensions, we expect that 

AF = 27F, (35) 

with 7 the line tension (the factor of two stems from the 
use of periodic boundary conditions, which lead to the 
formation of two interfaces in the system, see also the 
snapshot of Fig. \T7\ . Shown in Fig. [TBT a) is AF as a 
function of L; the data are indeed well described by a 
straight line through the origin. From the slope of the 
line, we obtain 7 = 0.158 k^T/a. Of course, we may also 
read-off the barrier in In (p, /i) , shown in Fig. llGf b) . 
Again, a linear increase of AF is observed, and from the 
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slope of the line we obtain 7 = 0.151 fceT/a, which is 
very close to our previous estimate. 

In summary, and in agreement with our theoretical re- 
sults, we find that an exponent p = 75 is high enough 
to induce a first-order transition in the off-lattice liquid 
crystal model of Eq. ([3]). The scaling of the coexisting 
densities with system size are well described by what is 
expected for such a transition. The same also holds for 
the growth of the free energy barrier AF. An interesting 
and unexpected result, which certainly requires further 
elaboration, is the scaling of the nematic order parame- 
ter, see Fig. [T31 If the trend continues for L ^ 00, long- 
range nematic order in a two-dimensional liquid crystal 
would, after all, be possible. 



B. Results for different p 

We have also measured the coexistence densities, 
chemical potential, and line tension for different values 
of p, while keeping 1/ = and e = 2.5. Here, we did not 
perform a detailed finite-size scaling analysis. Instead, 
the data for p ^ 75 were obtained in a single simulation 
of a (reasonably large) system. In Fig. [H we show the 
density of the isotropic phase, and of the nematic phase, 
as a function of p. Note that the region in between the 
curves corresponds to phase coexistence. If one performs 
a simulation in this region, for example by keeping the 
density fixed, snapshots strikingly reveal the two-phase 
coexistence, see Fig. [iTl Shown in Fig. [5] is the coex- 
istence chemical potential /i* versus p, compared to the 
DFT result. Finally, in Fig. El^b), we show the line ten- 
sion versus p. 



V. DISCUSSION AND SUMMARY 

In summary, we have provided strong evidence that off- 
lattice liquid crystals in two dimensions can also undergo 
first-order phase transitions. To this end, the pair poten- 
tial of Eq. ^ was introduced, constructed to be purely 
repulsive, short-ranged, and to obey inversion symmetry. 
The first-order transition takes place when the pair po- 
tential becomes sufficiently "sharp and narrow", i.e. for 
a sufficiently large exponent p in Eq. ([3]). A simple DFT 
calculation puts the threshold value at around p = 8. 
When p > 8, the theory predicts a first-order transition 
from a low-density (isotropic) phase, to a high-density 
(nematic) phase. In other words, there is a finite density 
gap between the two phases, as well as a jump in the 
nematic order parameter. 

An important conclusion of this work is that many of 
the trends predicted by the theory, also appear in the 
computer simulations. In other words, key properties of 
Eq. ([3]) are already captured at the mean-field level. This 
finding is somewhat surprising because, in low spatial 
dimension, mean-field is typically assumed to be unre- 
liable. The simulations, in agreement with the theory. 



find strong evidence of a first-order transition, already at 
p > 50. Note that we have backed our simulations with a 
detailed investigation of finite-size effects, and that these 
were shown to be consistent with a first-order phase tran- 
sition. We have not determined in our simulations the 
precise value of p where the crossover from continuous 
to first-order occurs, since such simulations would be ex- 
tremely time consuming, but we expect this will exceed 
the theoretical bound p — 8 significantly. 

Nevertheless, for those values of p where the simula- 
tions do observe the first-order transition, a profound 
density gap between the two coexisting phases is found, 
see the phase diagram of Fig. [4] In the limit of large p, 
theory and simulation are in qualitative agreement: both 
show an increase in density of the nematic phase with p, 
while the density of the isotropic phase is much less sen- 
sitive to p. At lower p, qualitative discrepancies arise, 
but these can be attributed to the large-p approximation 
used by the theory, which obviously breaks down here. 
Interestingly, when comparing the coexistence chemical 
potential /i* between theory and simulation, see Fig. [51 
the agreement is remarkably good, also at low p. Appar- 
ently, the breaking-down of the large-p approximation 
does not affect /i* as much as it does the coexistence 
densities. With regards to the fine tension, see Fig. [SI 
the agreement between theory and simulation is merely 
qualitative: by increasing p, the tension increases in both 
cases, but the actual numerical values diS^er profoundly. 
A possible explanation may be the presence of capillary- 
wave interface fluctuations. These fluctuations are ne- 
glected in the theory, while the simulation snapshot of 
Fig. [17] suggests that interface fluctuations are actually 
quite strong. 

A rather controversial result of this work is the value 
of the nematic order parameter S in the nematic phase. 
The DFT result, see Fig.[3l predicts rather large values of 
S", once the transition has become first-order. Of course, 
in realistic systems, one always has defects, which are ex- 
pected to destroy nematic order in the thermodynamic 
limit. Since such defects are discarded by the theory, we 
expect that Fig. [3] is merely an artifact, and that realis- 
tic systems in the thermodynamic limit will always have 
S — 0, regardless the value of p. Still, it is somewhat sur- 
prising that our computer simulations also suggest that 
5 > 0, see Fig. [T31 and that the finite-size effects in our 
data are not compatible with a decay of S to zero. At this 
point, the most likely explanation is that the decay of S 
with system size only shows-up in macroscopically large 
samples, which are clearly out of reach in any foreseeable 
simulation. 

Note that our theory has also made a number of 
intriguing predictions for the case where the transla- 
tional and orientational degrees of freedom are coupled, 
i.e. when v ^ 0. In this case, strong anchoring effects are 
predicted, as well as the possibility of the isotropic-to- 
nematic transition being pre-empted by freezing. Since 
our particles are ultrasoft, the formation of stable ag- 
gregate or cluster mesophases as encountered in various 
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FIG. 17: Simulation snapshot obtained in the coexistence region for p = 60. Each hne segment represents a particle. Clearly 
visible is that the system has phase-separated into an isotropic domain, and a nematic domain (where we leave open the question 
whether the nematic phase exhibits true long-range order, or only quasi long-range order). Note also that the interfaces are 
not flat, and that they appear to be decorated with capillary waves. 



soft-sphere systems |41|,|4^ is also possible. All these sce- 
narios need to be verified by computer simulation. We 
are currently developing new simulation methodology to 
study anchoring effects at the isotropic-nematic interface; 
some preliminary information about the method is al- 
ready available The application of these new tech- 
niques to the 2D liquid crystal model of the present work 
is therefore postponed to a future publication. 

Of course, it would be interesting if some of our find- 
ings could be confirmed in experiments. As we had 
already remarked in the Introduction, the condition of 
phase coexistence is obtained straightforwardly by keep- 
ing the density fixed at some value in the coexistence 
region. The problem will most likely be to achieve suffi- 
ciently "sharp and narrow" interactions. One possibility 
that we envision is to use a mixture of colloidal rods and 
non-adsorbing polymers. Despite the fact that colloidal 
rods cannot overlap, unlike the particles considered here, 
and that the addition of polymer renders the effective 
rod interactions attractive^ which may induce an addi- 
tional gas-liquid phase separation, we believe a first-order 
isotropic-nematic phase transition could be feasible. In 
particular, looking at thephase diagrams reported for 3D 
rod-polymer mixtures [44] we anticipate that the strong 
polymer-induced widening of the isotropic-nematic coex- 
istence region carries over to 2D systems as well, and 
might induce a first-order transition. A prerequisite for 
this scenario is that both the size ratio (of rod length to 



polymer radius of gyration) and the polymer concentra- 
tion are sufficiently large. 

Finally, we would like to remind the reader that the 
original idea of this work, namely to use "sharp and nar- 
row" interactions is not new, and goes back to the work of 
Ref. |l^. Here, the approximate correspondence between 
generalized XY-models and the Potts model [45| was al- 
ready exploited to demonstrate the possibility of having 
a first-order transition in a 2D spin system with contin- 
uous degrees of freedom. The extension of this work has 
been to apply the same ideas to off-lattice liquid crystals. 
Nevertheless, the approximate link to the Potts model 
can still be uncovered, using q cx p^l"^ [l^. Here, p is the 
exponent in Eq. and q the number of Potts states. 
Note that this relation is valid only asymptotically for 
large p. Many of the results of this work, for example 
the cross-over from a continuous to a first-order transi- 
tion (Fig. [1]), the variation of pL* with p (Fig. [5]), and 
even the increase of the line tension with p (Fig. [6]) , have 
their analogues in the 2D Potts model (see for example 
Ref. H^, where an explicit formula for the line tension of 
the Potts model is given). 
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